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CONTRACTOR REPORT 


THE M-INTEGRAL FOR COMPUTING STRESS INTENSITY FACTORS 
IN GENERALLY ANISOTROPIC MATERIALS 

1. INTRODUCTION 


The objective of this project is to develop and demonstrate a capability for computing stress 
intensity factors in generally anisotropic materials. These objectives have been met. 

The primary deliverable of this project is this report and the information it contains. In addition, 
we have delivered the source code for a subroutine that will compute stress intensity factors for 
anisotropic materials encoded in both the C and Python programming languages and made 
available a version of the FRANC3D program that incorporates this subroutine. 

Single crystal super alloys are commonly used for components in the hot sections of 
contemporary jet and rocket engines. Because these components have a unifonn atomic lattice 
orientation throughout, they exhibit anisotropic material behavior. This means that stress 
intensity solutions developed for isotropic materials are not appropriate for the analysis of crack 
growth in these materials. Until now, a general numerical technique did not exist for computing 
stress intensity factors of cracks in anisotropic materials and cubic materials in particular. Such a 
capability was developed during the project and is described and demonstrated herein. 

1.1 Background 

Expressions for the stress and displacement fields near the tip of a crack in an anisotropic 
material were first published by Sih, Paris, and Irwin [1965]. Their results were used by a 
number of authors to create numerical techniques for computing stress intensity factors in such 
materials (e.g., Atluri et al. [1975], Boone et al. [1987], and Tohgo et al. [1993]). However, a 
restriction not strongly emphasized in the original or any of the subsequent papers is that these 
results are only valid for the case where the plane of analysis (typically a plane perpendicular to 
the crack front) is a plane of material symmetry. 

The expressions for the stress and displacement field near a crack front in a generally anisotropic 
material were first published by Hoenig [1982]. These results do not appear to be widely 
referenced in the literature and it does not appear that a numerical technique for computing stress 
intensity factors based on these results has been developed previously. The Hoenig paper 
contains a significant number of typographical errors. We believe that these have all been 
identified and corrected in the work reported herein. 



1.2 Anisotropic Material Properties 


The primary materials of concern for the present study are single-crystal alloys. These materials 
exhibit cubic behavior, which means that they can be characterized by three properties (usually 
E, G, and v) and their orientation is given by three angles. The current implementation supports 
the more general case of arbitrarily oriented orthotropic materials, which are characterized by 
nine properties (E\, E 2 , £3, G 12 , G 23 , G 13 , V\ 2 , V23, and V] 3) and three angles. The fully general 
case of a material with 2 1 independent material properties could be accommodated easily in the 
current software, but as it is rare to assume such general behavior in engineering practice, this is 
not currently supported in order to simplify the interface. 

As mentioned, three angles are required to specify the orientation of the material axes relative to 
the Cartesian axes. There are a number of conventions for this. Here we adopted the use of Euler 
angles, which are commonly used in many finite element applications. Our convention conforms 
to that used by the ANSYS software, where the first rotation, 0 Z , is about the original z-axis. The 
second rotation, Q x , is about the rotated x-axis, and the third rotation, 9 V , is about the rotated y- 
axis. In all cases, a "right-handed" rule is used to determine the sign of the rotation. 

The remainder of this report is divided into seven sections. In the following section, the 
expressions for the stress and displacement fields near crack fronts in anisotropic materials are 
developed. In Section 3, the M-integral is fonnulated for anisotropic materials. In Section 4, the 
details of the present implementation of the M-integral are described. In Section 5, the results of 
verification studies with prescribed nodal displacements are presented. In Section 6, the results 
of verification studies of cracks in thick center-crack plates are presented. In Section 7, the 
results of a large number of analyses of cracks in single-crystal Brazilian disk specimens are 
summarized. Finally, Section 8 summarizes the report. 
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2. CRACK-FRONT STRESS AND DISPLACEMENT FIELDS 


2.1 General Anisotropy 

The stress and displacement fields in the vicinity of a crack-front in a general anisotropic 
material were presented by Hoenig [1982] based on the Lekhnitskii [1963] formalism. These are 
given by 






( 1 ) 


and 




( 2 ) 


where Re represents the real part of the quantity in parenthesis, two repeated indices obey the 
summation convention from 1 to 3, the coordinates x, y, z, r, and Prefer to the crack coordinates 
illustrated in Figure 1, and K/ are the stress intensity factors Ki, Kn, and Km- It is assumed that the 
stresses and displacements are functions only of the in-plane coordinates, r and 6. 



Figure 1. Local crack front coordinate system used throughout this report. 
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The quantities p t , A t , my, Ny, and Q t , are derived from the material compliance tensor, S ljk[ (where 
£y = Syki&ki)- The overbar indicates quantities specified in the Cartesian coordinate system. 

Before extracting the quantities, it is necessary to rotate the tensor to the crack-front coordinate 
system, Syu, and contract it to a 6x6 matrix, Sy, where i = 1..6. It is common in many branches of 
solid mechanics to work with the contracted form of the constitutive matrix. However, there are 
two different conventions for the ordering of the elements in the matrix. Both are used in this 
work. In the theoretical Sections (2 and 3), indices of the contracted matrix are taken in the order 
{x, y, z, yz, xz, xy} . This is consistent with the Hoenig and Sih papers and much of the literature 
for ortho tropic materials (principally composites). However, in the actual code and those 
portions of the report concerning the implementation, the contracted matrix indices are taken in 
the order {x, y, z, xy, yz, zx \ . This is consistent with most finite element software and much of the 
finite element literature. 

For plane stress conditions, the full compliance matrix is employed. For generalized plane strain 
conditions (which are assumed in the implementation and throughout this report), a reduced 
compliance matrix is obtained from 


, 6’ S’ 


(3) 


J 33 


where i,j = 1..6, S'y is symmetric, and S' a = S'y = 0. 

The parameters ju h i = 1..3, are the eigenvalues of the compatibility equations with positive 
imaginary parts. These are found from the characteristic sixth order polynomial equation 


l 4 {ju)l 2 {ju)-l:{ju)={) 


(4) 


where 

1,(m)=S’ ssM 2 - 2S' tsM + S' 4l 

iM=s;p-K+s' t y +(s' ls +s'J M -s'^ ( 5 ) 

iM=s;y- 2 sip + (2 s; 2 + 3 ; y- - 2 sy + s', , 

Solution of eqn. 4 leads to three pairs of complex conjugate roots. The three with positive 
imaginary parts are chosen. The parameters 2, are given by 

A i =- l Jf±\,i=\.3. ( 6 ) 

^2 \Mi / 


The parameters niy are given by 
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(7) 


>»„ = S'„ M ; -s' uM ,+s' n +z,(s; iMi -s; 4 ) 

m 2, = -S'y,+S' r Jfl +i(SL -S’,J /{/, ). 

m 3l = S'„ M , -S ' 6 +S: 2 /m, + Ate -S'JM,) 


The matrix AT is defined as 


N ij 


1 1 1 

-Mi ~Mi ~Mi 
- Aj — A 2 — Aj 


which has an inverse 


AT 1 = — 
' |7V| 


MiM mAi ^3 ^i Mi Mi 

JLl 3^1 /^i^3 >^3 ^i/3 ^i/j 

Z^i -^2 — Mi^i ^2 ~ Mi — Mi 


and |A] is the detenninate of Ny. Finally, 

Q t = a /cos 6* + //, sin# 


( 8 ) 


( 9 ) 


( 10 ) 


2.2 Special Anisotropies 

There are special anisotropies for which the above development is mathematically degenerate. 
These include monoclinic and cubic materials oriented such that z = 0 is a plane of material 
symmetry. Hoenig [1982] gives expressions for the stresses and displacements for the degenerate 
cases. However, the implementation follows the more widely referenced expressions given 
earlier by Sih et. al. [1965]. These stress and displacement expressions are 


<T= 
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Kj Re 
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Mi 
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and 



' 1 ] 

v Ml — Ml j 
' 1 ^ 

v Ml — Hi J 


( PiQi-PxQx ) 
( H2Q2 ~~ tf\Q\) 


( 12 ) 


In these expressions, ju\ and fii are the roots of the equation 

s;y - 2s;y + (2 s; 2 + s; ( y - 2 sy + s ’ 22 = 0 

with positive imaginary parts. 
jUi is the root of the equation 

S' 44 p 2 +2S' 45 p + S' 55 =0. 


with a positive imaginary part. 

The parameters pj and qj are defined as 

Pj =s^-s; 6 p j+ s: 2 , = i2 
Qj — S 12 /J ~ ^26 ^22 j/^ j 


Finally, 


Q : = yj cos^ + z/, sin# . 


(13) 


(14) 


(15) 


(16) 
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3. M-INTEGRAL FORMULATION 


3.1 The Relationship Between J and K 

As a first step in formulating an M - integral for anisotropic materials, the relationship between the 
./-integral and the stress intensity factors must be established. 

The J-integral developed by Rice [1968] is given by 


J = \[wn, 


-r 


dU: 


dx 


ds 


(17) 


where, for 2-D, the integration path Ais illustrated in Figure 2, W is the strain energy density 

(18) 


W = 2 (7 ij £ i j 


Hi is the outer nonnal to A A = u i is the displacement vector, and ds is the differential arc 
length along A 



Figure 2. The integration contour for the 2-D J-integral. 

By means of the crack closure integral, Hoenig [1982] developed the relationship between the J- 
integral and the stress intensity factors as 

J = V- [ A+K„ (!») 

where Im is the imaginary part of the bracketed quantities and summation is applied to repeated 
indices. 

For monoclinic or cubic materials oriented such that z = 0 is a symmetry plane, the relationship 
between J and K is 
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J = \ Im| 


s; 6 - 


M1+M2 


V MxMi J 


S’ 


22 


r i 4 

i S 22 


K,K n 

{ J 



+iJs'„su-s’iK 


(20) 


r 2 z^2 
III 


3.2 The M-Integral Expressions 

Yau et. al. [1980] developed the M-integral from the ./-integral as a way to extract the stress 
intensity factors for the three fracture modes from the global energy release rate. To obtain the 
M-integral from eqn. 19, two solutions are assumed and superposed. This is possible if the 
material is assumed to be linear elastic. Define 


cry = erf + ct f 

£j = 4" + e T 

(1) . (2) 
u , - u) 4- w. 


K. = K m + K 


( 2 ) 


( 21 ) 


Substitution of the stress, strain, and displacement fields into the J-integral expression (eqn. 17) 
leads to 


where 


J = J (1) +J (2) +M (1 ’ 2) 


( 22 ) 


J (]> =| 


du 

dx, 


(i) 


Ids 


i ) 


= f 1 

w (2 \ T t ; (2)5m,(> ] 

ds 


Jr( 

v J 



m (1 ’ 2) = 

f { W (h 2 ) n , T a) 8u ‘ } 

T d) 5 k, (1)> | 


1 ' dx 

1 

dx x J 


ds 


(23) 


The interaction strain energy density is given by 

W (1 ’ 2) =a m £ (2) =(T i2) £ m 

j v u y 


(24) 


and M* 12) is the M-integral. 

The superposed stress intensity factors in eqn. 21 can be substituted into eqn. 19 to yield 
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( 25 ) 


J = J m + J m Im K,JV:'/:f>) + A-; 21 Im(m 2 ,iV,:'^ [ ’) 

+ Im(m„JV,: l A:( 2) ) + Kf lm(m u N;lKf) 

+ k2> ) + K™ ’)] 


where 


J (1, =- 

J (2) = . 


#?> Im (m lt N:;Kf) + K$ Im (m v NjKf) + K>» V-'AT'”)] 

+AT!, 2) Im( mi ,JV- | A:f) + AT® Im^A^f )] 


. (26) 


Equating eqn. 22 and eqn. 25, and noting eqn. 23, leads to the definition of the M - integral as 

M =-i[<>im(». 2 ,iv: l ^ 2 >) + K™lm[m ll N:;Kf) 

+ K)» ) + K? Im(m u N:'K<>> ) . (27) 

+ K% Im (m„N-'Kf) +K% Im( m „ JV,; 1 Af < 2 > )] 


To this point, the J- and M-integrals have been presented as contour integrals. It is well known 
that in the context of finite elements the evaluation of area and volume integrals lead to more 
accurate and stable results than direct integration along a contour. Li et. al. (1985) introduced the 
2-D equivalent domain integral (EDI) form of the ./-integral that transforms eqn. 17 to an area 
integral. An analogous technique can be used to formulate a 3-D equivalent domain version of 
the M-Integral where the integration takes place over a volume. The resulting expressions are 


J (k) = f 

Jr 


-(h) 


du 


(k) 


dx. 


- W (k) 8 , 


v 


^-dV i = 1,2,3 k,j = 1,2 

jdXj 


M 


( 1 , 2 ) 


-L 


VI <T,/ 


(D ^ U i jg^(2) 


dx j 


- + cr,) 

,J dx , 


i j 


-!-dV i = 1,2,3 7 = 1,2 

jdXj 


(28) 


A C =J L q,ds. 


where 5] is the Kronecker delta and q is a function that is 1 at the stress intensity evaluation point 
and zero on the outer boundary of the domain of integration; q, is the value of the q function 
along the crack front. The domain of integration is a cylinder that encloses a portion of the crack 
front. This is illustrated in Figure 3. The q function can be thought of as a virtual crack 
extension. This is illustrated in Figure 4. 
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In order to evaluation the M-integral numerically, two solutions are chosen. Solution (1) is the 
finite element solution for the problem of interest. Solution (2) is an auxiliary solution; three 
solutions are chosen in order to obtain Kp\ Kjj l \ and AT///' 1 . These are denoted (2a), (2b), and 
(2c) and are chosen as follows 


Solution 

K, 

K n 

Km 

(2a) 

1.0 

0.0 

0.0 

(2b) 

0.0 

1.0 

0.0 

(2c) 

0.0 

0.0 

1.0 



Figure 4. The q function used in evaluating the 3-D M-integral. 


Substituting these into eqn. 27 yields the following expressions for M 


10 




M (h2a) = 
M"' U) = 
M' h2c> = 


2 K? Im (m 2l N :; ) + K'» )+ K'» + m 2i N :‘ )] 

K'," lm{m 2i N~l + m„N ) + 2 JC™ )+ K‘ V)' )] (29) 

Im(m 2l N;‘ + m 2i N2 ' ) + Kf Im(m„Ar,? + m,,^ 1 )+ 2 K% j 


and 


J 4 y a*! y a*, 2 


^-dA 


j dx i 


M {l ’ 2b) = JJ o-, (1) 


<3x, 


- + cr 


p M (i) 

(26) Jj /(l ’ 2b> S 

dx, lj 


dq 


J dx i 


dA 


m ( i ’ 2c> = r 

J 4 y cbc, y Sc, 12 


y dx j 


dA 


(30) 


These can be combined into one system of equations, 

2 \m(m 2i N -' ) Imfm,,^' 1 + m 2[ N ; 2 ) hn(m 3t N^ + m 2i N ^ ) 

lm(m 2i N; 2 + m u N-' ) 2 \m(m u N' i2 ) Im(m 3i 2V: 2 + m u AT 3 ) 

Im(m 2! 2V: 3 1 + m 3i N~' ) Im(m 1 ,JV /3 1 + m 3i N7 2 ) 2 Im(m 3; iV'' ) 


K<" 

K m 

K {1) 

^ in 


M a ' 2a) /A 
M^ ,lb) jA i 
m ( 1 ’ 2c) /4 


(31) 


that can be solved for the stress intensity factors. The terms in the coefficient matrix are obtained 
from eqs. 7 and 9 and depend on the material properties and orientation. 

Terms in the right-hand-side vector (eqn. 30) are obtained from the finite element results and 
eqs. 1 and 2. 

For the special case where z = 0 is a plane of material symmetry, the corresponding expressions 
for the M-integral are 


M 


M 


M 


(1,2a) 


= Im 


f 

S' x - 

f Mi + Mi 'l 

7 

S 22 



V 

l M \ M i ) 

y 

A 


M 1 M 2 ^ ii 


-S' 


MlMl ^ J 


K 


0 ) 


0,26) _ 


Im 


fj, x fj, 2 S n 


-S' 


22 

MiMi j 


«■? +((//, +m 2 )si< -s;,)k 


( 1 ) 

II 


(32) 


(1,2c) 


— / C f O ' _ O ' 2 V' (i) 

— y 3 44 ^ 55 ° 45 ^ III • 


These can be arranged into the system of equations 
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4 Ira 
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(1,2a) 


a; (1) 

,v // 
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4 

m (1 - 26) 

M (Uc) 


4. 


(33) 

where the terms in the right-hand-side vector are determined from eqn. 30 using eqs. 11 and 12 
to determine the stresses and displacements. 
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4. M-INTEGRAL IMPLEMENTATION 


4.1 Integration Domains and (/-Functions 

There are two important aspects of the implementation of the 3-D M - integral that do not follow 
directly from the formulation. These are the selection of the elements that will participate in the 
integration domain, and the related issue of selecting the virtual crack extension, the (/-function. 

Our experience with 2-D implementations of the ./-integral indicates that accurate stress intensity 
factors can be computed using eight triangular quarter-point elements around the crack tip. We 
obtain accurate results integrating only in these crack-tip elements. Other authors, however, 
report more accurate results when the integration is performed in "rings" of elements that are 
remote from, but completely encircle the crack tip (e.g., Banks-Sills and Sherman [1992]). In 3- 
D, unlike in 2-D, the domain of integration cannot be remote from the crack-front. That is, the 
crack-front elements must be included in the integration domain, if not, it cannot be 
demonstrated that the ./-integral (or M-integral) is path independent. 

Considering this experience and constraints, we have implemented the two options of evaluating 
the M-integral: either in the crack-front elements only (Fig. 5-8), or in the crack-front elements 
plus the first ring of elements surrounding the crack-front elements (Fig 9-12). The crack-front 
elements are assumed to be 15-noded wedge elements with the appropriate side nodes moved to 
the quarter-points. The elements in the next ring out from the crack front are assumed to be 20- 
noded brick elements. 

Related to the selection of elements in the integration domain is the selection of an appropriate 
virtual crack extension or (/-function. Banks-Sills and Sherman studied this issue in some detail 
[1988], They determined that the most accurate results are obtained if the integral is evaluated at 
the corner node of an element with a linear variation of the (/-function along the crack front in 
both sets of elements adjacent to the evaluation point. This (/-variation is illustrated in the right 
of Figures 5 and 9. For the purposes of this report, this is denoted as (/ (l) . This type of (/-function 
has the attractive feature that it can be used for nodes where the crack-front meets a free surface, 
as shown in Figures 6, 7, 10, and 11, which are denoted as cf~ ] and q {3) . 

An alternative (/-function that extends only the mid-side node is shown in Figures 8 and 12 
(denoted q (4} ). Banks-Sills and Sherman found that this extension was less accurate than the 
comer node extensions. However, this extension does overcome complications that arise for 
curved crack fronts. 

Curved crack fronts are usually approximated by piecewise linear segments. For anisotropic 
materials, one must consider the orientation of the material coordinate system relative to the 
crack front coordinate system. For cases q {2 \ c/ (3) , and (/ (4) the relative orientation of the material 
coordinate system can be determined from the tangent of the participating crack front segment 
and the local normal to the crack plane. For case q (X \ two crack front segments participate in the 
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integration domain, and, in general, their tangents will not be equal. This means that the local 
material coordinate system orientation is ambiguous. How this ambiguity is resolved in the 
current implementation is presented in Section 4.2.2. However, the accuracy of this resolution 
has not been assessed. 



Figure 5. The domain of integration and the virtual crack extension for q il) 

using one ring of elements. 



Figure 6. The domain of integration and the virtual crack extension for cf ] 

using one ring of elements. 
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Figure 7. The domain of integration and the virtual crack extension for q 

using one ring of elements. 




Figure 8. The domain of integration and the virtual crack extension for q 

using one ring of elements. 





Figure 9. The domain of integration and the virtual crack extension for q 

using two rings of elements. 



T z 

Figure 10. The domain of integration and the virtual crack extension for q {2) 

using two rings of elements. 



Figure 11. The domain of integration and the virtual crack extension for q (2) 

using two rings of elements. 



Figure 12. The domain of integration and the virtual crack extension for q (2) 

using two rings of elements. 
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A study of the relative accuracy of the combinations of element domains and ^-functions has 
been performed. The results are reported in Section 5. 

4.2 Detailed Implementation Steps 

The following sub-sections detail the steps in the current implementation of the M - integral for 
anisotropic materials. 

4.2.1 Compute material constitutive matrices 

Material properties are specified by nine material properties: E\, E 2 , £3, vn, v 2 3 , vn, Gn, G 23 , 
and G31. The orientation of the material coordinate system relative to the Cartesian coordinate 
system of the model is specified by three Euler angles, G z , d x , and 6 y . 

The material compliance matrix for an orthotropic material in the local material coordinate 
system is 

l/E l -v l2 /E 2 —v, 3 / E 3 0 0 0 

\/E 2 - v 23 / E 3 0 0 0 

1/E 3 0 0 0 

1/G 12 0 0 

1 / g 23 0 

sym 1 / G 23 

This is inverted to obtain the local material stiffness matrix 

C = S l . (35) 

A 3x3 transfonnation matrix that will rotate the local material coordinate axis to the material 
coordinate axes in the Cartesian space is detennined from the Euler angles as 

cos(YY ) cos(O v ) - sin(A ) sin(C. ) sin(A, ) - sin(A ) cos(A ) 

[/■] = sin(G_ ) cos(G, ) + cos(A)sin(Yf )sin(6 , v ) cos(C)cos(C ) 

-cos(/f)sin(C) sm(C) . (36) 

cos (A ) sin(A ) + sin( 0 . ) sin(Yf ) cos(Yf ) 
sin(YY ) cos (0 r ) - cos( 0 . ) sin(tY. ) cos(YY, ) 
cos(6> r )sin(6» v ) 

Components of the 3x3 rotation matrices are used to formulate 6x6 matrices that rotate the 
material stiffness and compliance matrices. These two matrices are different, and are given by 
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and 




f 2 

*11 

f 2 

*12 

* 1 . 3 

*i { t X2 

^ 12^13 

*13*11 

t 2 
*21 

t 2 
*22 

4 

*21*22 

*2.2*23 

*23*21 

1 2 
*31 

t 2 
*32 

t 2 
*33 

*31*32 

*32*33 

*33*31 

Of f 
^* 1 1*21 

Of t 
^“* 12*22 

Of t 

^* 13*23 

*11*22 + t l2 t 2 j 

*12*23 + *13*22 

*13*21 + *C 

Of f 
^* 21*31 

2*22*32 

Of t 

^* 23*33 

t 2X t 22 + t 22 t n 

*22*33 + *23*32 

*23*31 + *2V 

Of t 

. * 31*11 

Of t 

^* 32*12 

Of t 

^* 33*13 

*31*12 + *32*1 1 

*32*13 + *33*12 

*33*11 ^ 31 j 


4 

f 2 

l 12 

t 2 

l 13 

Of t 
^* 1 1*12 

Of t 

^* 12*13 

Of t 

^* 13*11 

*21 

1 2 

‘22 

1 2 

‘23 

Of t 
^* 21*22 

Of t 

^* 22*23 

Of t 

^* 23*21 

t 2 

*31 

f 2 

‘32 

f 2 

‘33 

Of t 

^* 31*32 

2*32*33 

Of t 

^* 33*31 

*11*21 

* 12*22 

*13*23 

*11*22 + *12*21 

*12*23 + *13*22 

*13*21 +*11*23 

*21*31 

*22*32 

*23*33 

*21*32 +*22*31 

*22*33 + *23*32 

*23*31 +*21*33 

*31*11 

*32*12 

*33*13 

*31*12 +*32*11 

*32*13 *33*12 

*33*11 +*31*13 


( 37 ) 


(38) 


The material matrices are then rotated into the Cartesian system by 

[c] = fc]|ck] r and [s]=fclsk] r - (39) 


4.2.2 Compute the crack-front coordinate system 

For cases q {2 \ q c '\ and q {4) , the local z-axis is equal to the tangent to the crack-front segment that 
is participating in the integration. The local y-axis is parallel to the local normal to the crack 
surface. The local x-axis is perpendicular to the local y- and z-axes. 

As mentioned above, the local crack front coordinate system for q {l) is ambiguous. In the current 
implementation, the local crack front axes are detennined from the coordinates of the nodes that 
are shared by the crack-front elements on either side of the evaluation location. This is illustrated 
in Figure 13. 



Figure 13. The definition of the local crack- front coordinates for case q {l \ based 

on specified nodal coordinates 
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4.2.3 Transform nodal coordinates, displacements, and material matrices 

The M-integral implementation requires that all input nodal coordinates, finite element 
displacements, and the material stiffness and compliance matrices be rotated to the local crack- 
front coordinate system. Furthermore, it requires that the nodal coordinates be translated so that 
the M-evaluation point is at the local coordinate origin. 

If the components of the normalized crack-front axes specified in the model Cartesian coordinate 
system are given by x\ y\ and z', then the 3x3 matrix that will transform nodal coordinates and 
displacements in the Cartesian system into the crack-front system is given by 


M- 


y, 

x' y' 

y X y y 

K y'y Z [ 


z x 

z[. 


(40) 


The corresponding transformations are 

{**} = [t']{x-x origin } and \u }= [t']{u}. (41) 

The material rotation matrices are similar to those given in eqs. 37 and 38, but use components 
of the t' matrix. The corresponding transformations are 

[C]=fi][c]k7 and [Skkpkf- < 42 > 


4.2.4 Determine the nature of the isotropy/anisotropy 

Components of the compliance matrix are examined to determine if the material is 1) isotropic, 
2) the local z = 0 plane is a plane of material symmetry, or 3) the material is generally 
anisotropic. 

For the isotropic and z = 0 symmetry cases, the compliance matrix will have distinct patterns of 
zero and non-zero entries, z = 0 symmetry will have the pattern 


[S] 


* * 
* * 
* * 
* * 

0 0 
0 0 


* * 
* * 
* * 
* * 

0 0 
0 0 


0 0 
0 0 
0 0 
0 0 

* * 
* * 


(43) 


The first check is to see if any of the upper right 6 entries of the compliance matrix are non-zero. 
If so, the material is assumed to be generally anisotropic. 
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An isotropic material has the pattern 



* 

* 

* 

0 

0 

0” 



* 

* 

* 

0 

0 

0 


[s] = 

* 

0 

* 

0 

* 

0 

0 

0 

0 

(44) 

* 

0 

0 


0 

0 

0 

0 

* 

0 



0 

0 

0 

0 

0 

* 


The second check is to see if the 14, 24, 

34, 

or 

56 elements are zero. If so, and the following 


conditions hold 

Su=S 2 2 = S 33 , S l2 = Sn = S 2 3 , S 44 = S 55 = S 66 , and S 44 = Sn/2(1« 0 (45) 

the material is assumed to be isotropic. Otherwise, z = 0 symmetry is assumed. 


4.2.5 Compute anisotropic parameters 

If the material has been determined to be z = 0 symmetric, then the p's are detennined by finding 
the roots to eqs. 13 and 14. 

If the material has been determined to be generally anisotropic, then the p's, A's, m's, and A r 1 are 
determined from eqs. 4 through 9. 


4.2.6 Compute the auxiliary nodal displacements 

For all elements in the domain of integration, we detennine the nodal displacements associated 
with the three auxiliary solutions, u {2a \ u (lh) , w l2cl . The generalized plane strain expressions for 
these values are given by eqs. 2 and 12 for the general anisotropy case and the z = 0 symmetry 
case, respectively. For isotropic materials in plane strain, the expression is 


u 

x 



1 

G\2n 


6 

cos— 

2 


l-2v , + sin 2 — 


V 


■J 


. e 

sin — 
2 


2 - 2v + cos 2 — 


0(~ „ 


sin — 
2 


2-2v-cos — 


V 


V 

e r 


■j 


■j 


cos— 

2 


e 


-l + 2v , + sin“ — 


V 


- J 


0 

0 

2 sin — 
2 


Kj 

K u 
I K„ 


.(46) 


These expressions give the displacements as a function of the in-plane coordinates r and 0. A 
general approach for determining these must account for curved crack fronts where the "ends" of 
the element may not be perpendicular to the crack front. The technique for doing this in the 
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current implementation is illustrated in Figure 14. The evaluation point x has element parametric 
coordinates //, and £ The associated point on the crack plane, x\ is determined from the 
parametric coordinates E, and 77 in the quadrilateral region that is the intersection of the 
evaluation ring with the crack plane. The associated point on the crack front, x", is determined 
from the parametric coordinate E, and the associated crack-front segment, r is the magnitude of 
the vector from x" to x, and 9 is the angle between this vector and the vector from x" to x'. 



Figure 14. A schematic of how the coordinates r and 9 are determined. 


4.2 . 7 Compute displacement derivatives 


du {l) du (2a) du (2b) du {2c) 

For each integration point within an element, we compute , , , and . 

dx dx dx dx 

This is done using standard finite element interpolation. The expression for the displacement 
derivatives is 


[du/dx] = [j ] _1 [dN/d£]{u } (47) 

where J is the Jacobian matrix and N represents the element shape functions. The Jacobian is 
computed from the product of the shape function derivatives with respect to the parametric 
coordinates and the element's nodal coordinates: 

[j] = [dN/d^]{x}. (48) 


4.2.8 Compute local stresses 

For each integration point, we compute the local element stresses. For the finite element stresses, 
the expression is 
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( 49 ) 


°\v 


du x /dx 



du y fdy 

a xy 

• = [4 

du, / dz 

du x /dy + du y /dx 

CT VZ 


du y /dz + dujdy 



du,/dx + du x !dz 


For the auxiliary solutions (c/ ?a) , o* 2b> , and c/ 7c) ), stresses are given by eqn. 1 for the general 
anisotropy case, eqn. 11 for z = 0 symmetry, and for the case of isotropy by 


<j. 


cr 


a 


>’ z 


cr. 


V2 


m 


9 

cos— 

2 


V 

e r 


cos— 

2 


' . 9 . 3 9 

1 - sin— sin — 
2 2 

, . 9 . 39 

1 + sin— sin — 
2 3 


V ^ J J 

. 9 9 39 

sin— cos— cos — 

2 2 2 

0 

0 


9 


f 


-sin — 
2 


9 29 


2 + cos— cos — 

V 2 2 y 

. 0 9 39 

sin— cos— cos — 

2 2 2 


9 

cos— 

2 


f 


, . 9 . 39 

1 - sin — sin — 
2 2 

0 

0 


0 

0 

0 

9 

cos— 

2 

. 0 
-sin — 
2 


Kj 

K„ 


K 


m 


.(50) 


For all cases, cr~ comes from the plane strain condition 

= (S 3l <r x + S 32 a y + S, 4 a xy + S 35 cr yz + S 36 a zx )/S 33 . (5 1) 

4.2.9 Compute the contribution to the M-integral 

The M-integral is evaluated numerically using Gauss integration. It has contributions from all 
Gauss points for all elements in the integration domain. The expressions for the numerical 
integration is 


ttelem #gp 

elem = 1 gp = 1 
#elem #gp 

M il2h) = £ t WgP \j ; 

elem=\ gp = 1 
#elem #gp 

= z 

elem = 1 gp = 1 


gP 


a ( 1 )^ +CT ( 2 «)^!_ ff ( 1 , 2 ^ 

dx 1 dx i 


du 


( 26 ) 


(7 ? * — — 
V dx \ 


- + <7 , 


( 26 ) 


du; 


0 ) 


dx. 


i W {l2h) S, 


dq 

jdxj 

'g- (52) 

J dX J 


r 


SP 




du { , 2c) ( 2 c) du (X) 


dx j 


■ + (7 ii 


dx j 


i 




where vv^, is the Gauss integration weight and \J gp \ is the determinate of the Jacobian matrix 
evaluated at the integration point. Summation convention is assumed with i= 1, 2, 3 and j = 1, 2, 3. 
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The current implementation uses 27 integration points for brick elements and 28 integration 
points for wedge (crack-front) elements. 


4.2.10 Compute the area of the virtual extension 

A q is the area of the virtual extension along the crack front. For cases c/ 2 \ c / (3) , and q {4) only one 
crack-front segment participates in the integration. If this is denoted /, and, noting that the virtual 
extension at the evaluation point is 1 , 

C=C=V and A T=V (53) 

For case q 1 \ two crack- front segments participate in the integration. If these are denoted l\ and 
h, then 


4.2.11 Compute the stress intensity factors 


(54) 


Finally, the stress intensity factors are computed. Eqn. 31 is used for general anisotropy, eqn 33 
for z = 0 symmetry, and for isotropy, the expression is 


i'-s) 

E 

0 

0 


0 

2 M) 

E 

0 


0 

0 

2(1 + v) 


Kf 


M {X ' 2a) / A q 

K u 

> = < 

M«’ 2b) /A q 

K ni , 


/ A q 


(55) 
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5. PRESCRIBED DISPLACEMENT TESTS 


An initial series of tests were performed to verify the implementation. In these tests, no finite 
element analysis was performed. Rather, nodal displacements for all elements participating in the 
M-integral evaluation were determined from eqs. 2, 12, and 46, assuming unit stress intensity 
values. The M-integral is evaluated, and the computed stress intensity values are compared to the 
prescribed values. 

The purpose of these tests is two-fold. First, as verification; if the M-integral evaluation is 
encoded correctly, then the computed stress intensity factor should be close to the prescribed 
value. Second, the absolute difference between the prescribed and computed values give a 
quantitative assessment of the "intrinsic" error in the computations and an indication of the 
relative performance of differing q functions and numbers of element rings. By intrinsic error, 
we mean those errors that arise due to the approximate numerical techniques used to evaluate the 
integrals. 

Table 1 shows the computed stress intensity factors for isotropy, z = 0 symmetry, and general 
anisotropy. For each, all three fracture modes are examined for both 1 and 2 element rings 
around the crack front for four different q functions. 

Of all the tests, the largest error was 4.7%. This was for general anisotropy, mode I, one ring of 
elements, and q {4) . In general, all the q i4] tests with one ring of elements perfonned poorly. If 
these are discounted, then the largest error for the one-ring tests was 0.89%. Among the two- 
ring tests the largest error was 0.46%. 

Based on these tests, the following conclusions are drawn: 

1. All combinations of element rings and q functions produce accurate results with the 
exception of q {4) with one ring of elements. This configuration should be avoided. 

2. In general the two-ring configurations are slightly more accurate than the one-ring 
configurations, but the difference is not significant. 
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6. THICK PLATE ANALYSES 


A series of finite element analyses were performed to verify the perfonnance of the M- integral 
calculations. The models were essentially very thick center-cracked plates. The objective was to 
make the plates thick enough that plane-strain conditions would be approximated near the center 
of the plate. 

The geometry of the model is shown in Figure 15. The plate is simply supported as shown in 
Figure 16a. Two load cases, mode I (Figure 16b) and mode II (Figure 16c) were analyzed. All 
length and force units are normalized with respect to those used for the material properties. 


Table 1. Computed stress intensity factors for prescribed nodal displacements. 





prescribed K : 



prescribed K u 


prescribed K m 

rings 

(n) 

<r 

K, 

K n 

Km 

K, 

K u 

K m 

K, 

K n 

K m 

Isotropic 


1 

1.0089 

0.0 

0.0 

0.0 

1.0062 

0.0 

0.0 

0.0 

1.0070 

1 

2 

1.0089 

0.0 

0.0 

0.0 

1.0062 

0.0 

0.0 

0.0 

1.0070 

3 

1.0089 

0.0 

0.0 

0.0 

1.0062 

0.0 

0.0 

0.0 

1.0070 


4 

0.9591 

0.0 

0.0 

0.0 

0.9647 

0.0 

0.0 

0.0 

0.9562 


1 

0.9984 

0.0 

0.0 

0.0 

0.9987 

0.0 

0.0 

0.0 

0.9992 

9 

2 

0.9984 

0.0 

0.0 

0.0 

0.9987 

0.0 

0.0 

0.0 

0.9992 

Z 

3 

0.9984 

0.0 

0.0 

0.0 

0.9987 

0.0 

0.0 

0.0 

0.9992 


4 

0.9984 

0.0 

0.0 

0.0 

0.9987 

0.0 

0.0 

0.0 

0.9992 

Z symmetry 


1 

1.0080 

0.0 

0.0 

0.0 

1.0065 

0.0 

0.0 

0.0 

1.0069 

1 

2 

1.0080 

0.0 

0.0 

0.0 

1.0065 

0.0 

0.0 

0.0 

1.0069 

3 

1.0080 

0.0 

0.0 

0.0 

1.0065 

0.0 

0.0 

0.0 

1.0069 


4 

0.9534 

0.0 

0.0 

0.0 

0.9624 

0.0 

0.0 

0.0 

0.9543 


1 

1.0034 

0.0 

0.0 

0.0 

1.0009 

0.0 

0.0 

0.0 

1.0001 

9 

2 

1.0034 

0.0 

0.0 

0.0 

1.0009 

0.0 

0.0 

0.0 

1.0001 

z 

3 

1.0034 

0.0 

0.0 

0.0 

1.0009 

0.0 

0.0 

0.0 

1.0001 


4 

1.0034 

0.0 

0.0 

0.0 

1.0009 

0.0 

0.0 

0.0 

1.0001 

general anisotropy 


1 

1.0078 

0.0 

- 0.0007 

0.0 

1.0064 

0.0 

- 0.0013 

0.0 

1.0073 

1 

2 

1.0078 

0.0 

- 0.0007 

0.0 

1.0064 

0.0 

- 0.0013 

0.0 

1.0073 

3 

1.0078 

0.0 

- 0.0007 

0.0 

1.0064 

0.0 

- 0.0013 

0.0 

1.0073 


4 

0.9529 

0.0 

- 0.0024 

0.0 

0.9702 

0.0 

- 0.0062 

0.0 

0.9591 


1 

1.0046 

0.0 

0.0018 

0.0 

0.9968 

0.0 

0.0045 

0.0 

0.9980 

9 

2 

1.0046 

0.0 

0.0018 

0.0 

0.9968 

0.0 

0.0045 

0.0 

0.9980 

Z 

3 

1.0046 

0.0 

0.0018 

0.0 

0.9968 

0.0 

0.0045 

0.0 

0.9980 


4 

1.0046 

0.0 

0.0018 

0.0 

0.9968 

0.0 

0.0045 

0.0 

0.9980 
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Figure 15. The geometry used for the thick plate analysis. 



Figure 16. The simply supported boundary conditions (a), and the mode I (b) 
and II (c) load cases used in the analyses. 

Three sets of material properties were investigated 

Isotropic: 

£= 1000.0, v= 0.25 

Cubic material oriented with z = 0 symmetry: 

E\ = £2 = £3 = G \2 = C23 = C13 = 1000.0 
V12 = V23 = V13 = 0.25 
6 Z = 45.0°, G x = 6 y = 0.0° 
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Cubic material oriented with a general orientation: 

E\ = Ei — Et, = G \2 = G 23 = G \2 = 1000.0 
V\2 = V23 = Vl3 = 0.25 
ft = ft = ft = 45.0° 

The model and meshes were generated with FRANC3D. The finite element analysis was 
performed with ANSYS. The resulting nodal displacements were read back into FRANC3D 
where the M-integrals were evaluated to determine the stress intensity factors. 

The same mesh was used for all analyses; it is shown in Figure 17, with a detail of the crack 
region shown in Figure 18. The mesh is predominantly 10-noded tetrahedral elements. 15-noded 
quarter point wedge elements are used around the crack front. One ring of 20-noded brick 
elements is used around the crack-front elements. 13-noded pyramid elements are used to 
transition from the brick elements to the tetrahedral elements in the bulk of the mesh. 

Because of the simple loading and the simple supports, the well-known K solution for a center- 
cracked plate is valid for all three material behaviors. The relative errors between the theoretical 
solution and the numerical solutions are shown graphically in Figure 19. It is expected that there 
will be a relatively large error near the surfaces of the plate where plane strain conditions are not 
satisfied. It can be seen that near the center of the plate the M-integral solutions are very 
accurate. 



Figure 17. The mesh used for the thick plate verification analyses. 
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Figure 18. A detail of the thick plate mesh showing the crack region. 



LI 


-isotropic Mode I 

- z=0 sym Mode I 

- general Mode I 

- isotropic Mode II 
-z=0 sym Mode II 
-general Mode II 


normalized distance along crack front 


Figure 19. The relative errors between the theoretical and numerical solutions for the stress 
intensity factors for a thick center-cracked plate. 
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7. BRAZILIAN DISK ANALYSIS 


A large number of analyses were perfonned to detennine a stress intensity solution for a 
Brazilian disk test specimen. The specimen geometry is shown in Figure 20. The material 
properties considered are those for a typical single crystal nickel alloy used in an engine hot 
section. An isotropic material was considered as well. In testing, it is observed that the two ends 
of the crack do not usually grow at the same rate, therefore, a range of a\ and ai lengths are 
considered. 



Figure 20. The geometry of the Brazilian disk test specimen. 


The single crystal material has cubic symmetry. This means that E = fsooi = £010 = Tmo, v = 
M)oi/oio = Km/ioo = Hmo/ioo, and G = Gooi/oio = GW/ioo = Goio/ioo- The properties used here are E/G 
= 0.952 and v= 0.4. The properties used for the isotropic material are E/G = 2.8 and v= 0.4. 

Two different material orientations were considered for the anisotropic material, <11 1/1 J_0> and 
<111/11 2> where <y/x> defines the orientation of the local crack front axes relative to the 
material axes. These orientations are shown graphically in Figure 21, where the red axes indicate 
the local crack front-axes. The corresponding Euler angles are 9 Z = 50.7685, 9 X = -24.0948, and 
9 y = 26.5650, for the < 1 1 1/1 j_0> case, and 9 Z = -35.2644, 9 X = -45.0, and 9 V = -90.0, for the 
<11 1/1 12> case. 
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Figure 21. A graphical illustration of the cubic material orientation relative to the local crack- 

front axes. 

Twenty-eight crack configurations were considered. These are shown in Table 2. Eleven load 
orientations, Os, were considered. The isotropic and <111/1J_0> material configurations give 
symmetric results for negative and positive Wangles. For these cases, 0°, 4°, 8°, 12°, 16° and 20° 
were analyzed. For the <111/11 2> material orientation, -4°, -8°, -12°, -16°, and -20° were 
analyzed in addition to the positive angles. 

Table 2. Analyzed crack configuration. 


2a\lW 

2 a 2 /W 

0.2 

0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8 

0.3 

0.3, 0.4, 0.5, 0.6, 0.7, 0.8 

0.4 

0.4, 0.5, 0.6, 0.7, 0.8 

0.5 

0.5, 0.6, 0.7, 0.8 

0.6 

0.6, 0.7, 0.8 

0.7 

0.7, 0.8 

0.8 

0.8 


A typical mesh used in these analyses is shown in Figure 22 (with greatly magnified 
displacements). The meshes were two elements thick and periodic boundary conditions were 
applied to the front and back surfaces to enforce plane strain behavior. The meshes were 
predominantly 15-noded wedge elements, with quarter-point wedge elements used at the crack 
fronts and 20-noded brick elements used in the first ring of elements surrounding the crack-front 
elements. The M-integral was evaluated using q (l) . 
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Figure 22. A typical mesh used in the disk analyses (with greatly 
magnified displacements). 


Typical analysis results for are shown graphically in Figure 23. Complete analysis results are 
given in graphical and tabular form in the appendices. 
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Figure 23. Typical Brazilian disk analysis results, where 
K n nWt . , 

F n = — . ” , with n = I, II. 
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8. SUMMARY 


This report describes the successful effort to develop and verify a numerical technique for 

extracting stress intensity factors from finite element results for anisotropic materials. The report 

can be summarized as follows: 

• Expressions for the stress and displacement fields near a crack front in anisotropic material 
are presented. 

• The expressions for the M-integral that can be used to determine the mode-wise stress 
intensity factors in such materials are developed. 

• The details of an implementation of the anisotropic M-integral are presented. 

• The M-integral is shown to give consistent and accurate results for test cases where nodal 
displacements are prescribed to be the theoretical near crack-front values. 

• The M-integral is shown to give accurate results for a simple model with a known stress 
intensity factor solution. 

• The new capability is used to develop the /f-solution for crack growth in an anisotropic 
Brazilian disk specimen. 
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APPENDIX A: BRAZILIAN DISK ANALYSIS RESULTS 


This appendix contains stress intensity factor results for the Brazilian disk specimens. 

Analyses were perfonned for a cubic material with properties E/G = 0.952 and v= 0.4 and for an 
isotropic material with properties E/G = 2.8 and v= 0.4. Two material orientations relative to the 
crack-font were considered for the cubic material, <111/1J_0> and <111/11 2 > where <v/x> 
defines the orientation of the local crack front axes relative to the material axes. 

For the isotropic material and the cubic material in the <1 1 1 / 1 J_0> orientation, Wangles of 0°, 4°, 
8°, 12°, 16° and 20° were analyzed. Due to material symmetries, these results are symmetric 
about 0= 0°. The cubic material in the <111/1 12> orientation is not symmetric about 6= 0°. For 
that case, -4°, -8°, -12°, -16°, and -20° were analyzed in addition to the positive angles. 

A total of 644 analyses were performed. 

The first part of this appendix contains the stress intensity factor values in tabular form. Each 
page contains the data for one load orientation for one material type/orientation. 

The second part of this appendix contains plots of nonnalized stress intensity factors for Modes I 
and II. The normalized stress intensity factors are defined as 

K nWt 

F n = — , " where n = I, II. 

P^]27r(a 1 + a-,) 
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tropic materials using the M-integral. Included are examples of this solution applied to Brazilian 
disk crack growth specimens. 
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